home *** CD-ROM | disk | FTP | other *** search
- #include <stdio.h>
- #include <stdlib.h>
- #include <string.h>
- #include <ctype.h>
- #include "foil_lib.h"
- #include "lib.h"
-
- /*
- * get_coords
- */
- static int get_coords(FILE *fp, float *pt1, float *pt2, float *pt3)
- {
- char buf[TXTSIZ], *p;
- float x1, x2, x3;
- int n;
-
- while (1) {
- if(!fgets(buf, TXTSIZ, fp))
- return EOF;
-
- /* skip comments */
- p = buf;
- while(isspace(*p)) p++;
- if (in(*p, "%#\n"))
- continue;
-
- /* scan string */
- n = sscanf(buf, "%f %f %f", &x1, &x2, &x3);
-
- if (n == EOF)
- return EOF;
-
- /* valid data? */
- if (n > 0)
- break;
-
- /* otherwise who knows? just keep reading */
- }
-
- *pt1 = x1;
- *pt2 = x2;
- if (pt3)
- *pt3 = x3;
- return n;
- }
-
- int read_foil(FILE *fp, char *f, airfoil_data_t *info)
- {
- static float X[MAXPOINTS]; /* x points */
- static float Y[MAXPOINTS]; /* y points */
- static float YY[MAXPOINTS]; /* other y, for NACA format */
- float x=0., y=0., yy=0., minx, maxx, miny, maxy;
- int i, j, n, nread, npoints;
- char *p;
-
- fgets(info->title, TXTSIZ, fp);
- if((p = strchr(info->title, '\n')) != 0)
- *p = 0;
- n = get_coords(fp, &x, &y, &yy);
-
- if(n < 2 || n > 3) {
- fprintf(stderr, "Syntax error in line 2 of \"%s\"!\n", f);
- return 0;
- }
-
- minx = maxx = X[0] = x;
- Y[0] = y; YY[0] = yy;
- i = 1;
-
- while(n == (nread = get_coords(fp, &x, &y, &yy))) {
- X[i] = x; Y[i] = y; YY[i] = yy;
- minx = min(x, minx);
- maxx = max(x, maxx);
-
- if(i++ >= MAXPOINTS) {
- fprintf(stderr,
- "Too many points! Change MAXPOINTS (currently = %d) and recompile.\n", i);
- return 0;
- }
- }
-
- /* any error codes? could be more robust */
- switch (nread) {
- /* data OK */
- case 0: break;
- /* EOF */
- case -1: break;
- /* bad input */
- case -2: fprintf(stderr, "Syntax error at \"%s\" (line %d)\n", f, i+2);
- break;
- /* bogus - shouldn't get here */
- default: /* NOTREACHED */
- break;
- }
-
- if(n == 3) {
- /*
- * NACA format; combine into one array.
- */
- npoints = 0;
- for(j = i-1; j > npoints; j--, npoints++) {
- float temp;
- temp = X[npoints]; X[npoints] = X[j]; X[j] = temp;
- temp = Y[npoints]; Y[npoints] = Y[j]; Y[j] = temp;
- }
- npoints = i;
- for(j = 1; j < i; j++, npoints++) {
- X[npoints] = X[i];
- Y[npoints] = YY[i];
- }
- }
- else
- npoints = i;
-
- if(X[0] != X[npoints-1] || Y[0] != Y[npoints-1]) {
- fprintf(stderr,
- "Warning: data does not start and end at the same point.\n");
- X[npoints] = X[0]; Y[npoints] = Y[0];
- npoints++;
- }
- if(minx < 0.0)
- for(i = 0; i < npoints; i++)
- X[i] -= minx;
-
- if(maxx == 0.0) {
- fprintf(stderr, "Error: no x values?\n");
- return 0;
- }
-
- info->npoints = npoints;
-
- miny = maxy = 0.;
- for(i = 0; i < npoints; i++) {
- info->points[i].x = X[i]/maxx;
- info->points[i].y = Y[i]/maxx;
- maxy = max(maxy, info->points[i].y);
- miny = min(miny, info->points[i].y);
- }
- info->maxy = maxy;
- info->miny = miny;
-
- return 1;
- }
-
- void writefoil(FILE *fout, airfoil_data_t *foil, int tchanged, int cchanged)
- {
- int i;
-
- fprintf(fout, "%s", foil->title);
- if(tchanged)
- fprintf(fout, " T-%.3g%%", foil->thickness * 100);
- if(cchanged)
- fprintf(fout, " C-%.4g%%", foil->camber * 100);
- fputc('\n', fout);
-
- for(i = 0; i < foil->npoints; i++)
- fprintf(fout, "%f %f\n", foil->points[i].x, foil->points[i].y);
- }
-
- /*
- * Interpolate the function at x, given 2 points on either side. If
- * straightl is set, just do a straight line (for the first intervals, when P4
- * or P1 may not exist).
- *
- * If P1 and P4 are on opposite sides of (P2,P3) we use the straight line.
- */
- #define lin_interp(x, x2, y2, x3, y3) ((y2)+((x)-(x2))/((x3)-(x2))*((y3)-(y2)))
- static double interp(double x,
- double x1, double y1, double x2, double y2,
- double x3, double y3, double x4, double y4, int straightl)
- {
- double i1x, i1y, i2x, i2y; /* the two interior points for interpolation */
-
- /* A is (P1,P2); B is (P2,P3); C is (P3,P4). */
- point_t ma, mb, mc;
- ma.x = x2 - x1; ma.y = y2 - y1;
- mb.x = x3 - x2; mb.y = y3 - y2;
- mc.x = x4 - x3; mc.y = y4 - y3;
-
- if(straightl ||
- ((m_gt(ma, mb) && m_gt(mb, mc)) ||
- (m_gt(mc, mb) && m_gt(mb, ma))
- )) {
- if((x3 - x2) < EPSILON)
- /* Indeterminate */
- return y3;
- return lin_interp(x, x2, y2, x3, y3);
- }
-
- i1x = x2 + (x3-x2)/3.;
- i2x = x2 + (x3-x2)*2./3.;
- i1y = lin_interp(i1x, x1, y1, x2, y2);
- i2y = lin_interp(i2x, x3, y3, x4, y4);
- if(x < i1x)
- return lin_interp(x, x2, y2, i1x, i2y);
- if(x > i2x)
- return lin_interp(x, i2x, i2y, x3, y3);
- return lin_interp(x, i1x, i1y, i2x, i2y);
- }
-
- /*
- * Takes an array of points, and for each each one finds the corresponding
- * point on the other side by interpolating. Returns corresponding points in
- * the "new" array.
- *
- * For now, this cheap version just finds the linear interp on the other side.
- * Most airfoil data seems pretty smooth and the data points are close
- * together so the linear interpolation works fine; so let's leave the "right"
- * interpolation for another time. "TODO 3.?" I guess.
- */
- void get_otherside(point_t *new, point_t *old, int n)
- {
- int i, j, do_next, other;
- /* These are for the flags: si means i is at the start; sj is j at start. */
- int s, si=1, sj=1;
- int nnl, nnr, neighbour;
-
- /* We run the loop from both sides until they meet */
- new[0].x = old[0].x;
- new[0].y = old[0].y;
- new[n-1].x = old[n-1].x;
- new[n-1].y = old[n-1].y;
-
- i = 0; j = n-1;
-
- while(i < j) {
-
- if(old[i+1].x < old[j-1].x) {
- j--; do_next = j; other = i;
- neighbour = i+1; nnl = i-1; nnr = i+2;
- sj = 0; s = si;
- }
- else {
- i++; do_next = i; other = j;
- neighbour = j-1; nnl = j+1; nnr = j-2;
- si = 0; s = sj;
- }
- #ifdef CheapInterp
- new[do_next].x = old[do_next].x;
- new[do_next].y = lin_interp(new[do_next].x, old[other].x, old[other].y,
- old[neighbour].x, old[neighbour].y);
- #else
- new[do_next].x = old[do_next].x;
- new[do_next].y = interp(new[do_next].x,
- old[nnl].x, old[nnl].y,
- old[other].x, old[other].y,
- old[neighbour].x, old[neighbour].y,
- old[nnr].x, old[nnr].y, s);
- #endif
- }
-
- #ifdef Debug
- fputs("Interpolation debug\n", stdout);
- for(i = 0; i < n; i++)
- printf("%.5f %.5f %.5f %.5f\n", new[i].x, new[i].y,
- old[n-i-1].x, old[n-i-1].y);
- #endif
- }
-
- /*
- * Very similar to get_otherside(); instead of interpolating in the other
- * side, we interpolate in a different array.
- *
- * Caveat: when the curve turns the corner, we need to reverse the sense, so
- * we fill in the points in two parts: first the upper, then the lower.
- */
- void fill_in(point_t *interped, int n1, point_t *pts, int n2)
- {
- int i, j=0;
-
- for(i = 0; i < n1; i++) {
- if(pts[j+1].x > pts[j].x || interped[i+1].x > interped[i].x)
- break;
- while(pts[j+1].x > interped[i].x)
- j++;
- #ifdef CheapInterp
- interped[i].y = lin_interp(interped[i].x, pts[j].x, pts[j].y,
- pts[j+1].x, pts[j+1].y);
- #else
- /* ... *** not done yet *** ... */
- #endif
- }
- for(; i < n1; i++) {
- while(pts[j+1].x < interped[i].x)
- j++;
- #ifdef CheapInterp
- interped[i].y = lin_interp(interped[i].x, pts[j].x, pts[j].y,
- pts[j+1].x, pts[j+1].y);
- #else
- /* ... *** not done yet *** ... */
- #endif
- }
- }
-